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We study the effect of disorder on the propagation of collective excitations in a disordered Bose 
superffuid. We incorporate local density depletion induced by strong disorder at the meanfield level, 
and formulate the transport of the excitations in terms of a screened scattering problem. We show 
that the competition of disorder, screening, and density depletion induces a strongly non-monotonic 
energy dependence of the disorder parameter. In three dimensions, it results in a rich localization 
diagram with four different classes of mobility spectra, characterized by either no or up to three 
mobility edges. Implications on experiments with disordered ultracold atoms are discussed. 
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I. INTRODUCTION 

The dynamics of correlated quantum systems attracts 
a growing attention sparked by the recent development of 
quantum devices with long coherence times and dynam¬ 
ical control of parameters, e.g. superconducting circuits 
and ultracold atoms [1] . An additional asset of the latter 
is that disorder may be introduced in a controlled way [2] . 
Disorder may strongly affect dynamical processes, mainly 
due to Anderson localization [3]. Understanding the in¬ 
terplay of disorder and interactions in dynamical quan¬ 
tum systems is thus of fundamental importance and lo¬ 
calization in quantum systems is still the subject of active 
research [4-11]. This topic has also been addressed in the 
context of quasi-periodic systems [12, 13]. 

In correlated quantum systems, most basic dynami¬ 
cal processes are determined by the transport properties 
of their collective excitations [14]. An important start¬ 
ing point in the understanding of localization in corre¬ 
lated systems thus relies on classification according to 
the symmetries of their excitations [15]. For Fermi sys¬ 
tems, it is mostly based on the three classes of random 
matrices [16] as well as chiral or particle-hole symme¬ 
tries [17]. For Bose systems, a strong distinction arises 
between Goldstone and non-Golstone modes [18]. For in¬ 
stance, in a Bose superfluid, while localization is at its 
strongest at low energy for particle-like excitations, it is 
suppressed for phonon excitations [5, 7, 10]. This con¬ 
clusion is based on a weak disorder analysis and holds 
in dimension d < 2 where localization occurs for arbi¬ 
trary weak disorder. It is, however, challenged in higher 
dimension where the onset of the Anderson transition re¬ 
quires sufficiently strong disorder, which may alter the 
very nature of the excitations. 

In this work we study the transport of collective excita¬ 
tions in a disordered, weakly-interacting Bose superfluid 
in dimension higher than one. We show that the compe¬ 
tition of disorder, screening, and density depletion yields 
a strongly non-monotonic and non-universal energy de¬ 
pendence of the disorder parameter, which controls the 
localization properties. In three dimensions (3D), our 
analysis indicates that the localization diagram exhibits 



Figure 1. One-dimensional cut of the density profile of a Bose 
superfluid in a disordered potential: Exact numerical solution 
of the GPE (1) [shaded area, nc(r)] vs. self-consistent solution 
using Eq. (2) [dashed line], for a disordered potential [full line, 
U(r)] of amplitude Ur = 0.87/i and correlation length ctr = 
The density profile follows the modulations of the potential, 
smoothed at the length scale of G and may be locally depleted 
around disorder maxima. 


several classes of mobility spectra, characterized by either 
no or several mobility edges. We finally discuss implica¬ 
tions of these localization properties on quantum-quench 
experiments with disordered ultracold atoms. 


II. TRANSPORT THEORY OF COLLECTIVE 
EXCITATIONS 

A. Meanfleld scattering theory 

To study the transport of collective excitations in the 
presence of disorder, it is worth devising a scattering 
problem. For weakly-interacting Bose superfluids, we 
may rely on meanfield theory [19]. The background 
density field nc(r) obeys the Gross-Pitaevskii equation 
(GPE), 


[- fi^V'^/2m + V{r) - fi + gnc{r)]y/n^ = 0, (1) 

where m is the particle mass, /r is the chemical poten¬ 
tial, and g > 0 is the coupling constant of short-range 
repulsive interactions. The disordered potential U(r) is 
chosen to be spatially homogeneous, isotropic, and of 
vanishing statistical average. For weak disorder, Eq. (1) 
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can be solved perturbatively [20, 21]. Below, we con¬ 
sider regimes of strong disorder and it is necessary to 
extend this approach, including possible local depletion 
of the density around the disorder maxima. To do so, we 
generically write the density field in the form 

nc(r) = [/X-? 7 (r)-h AJ/g, (2) 

where the field 77 (r) describes the modulations of the den¬ 
sity due to the disorder, and the quantity A is a shift in 
the chemical potential. The latter allows us to impose the 
conventional condition that 77 (r) is of zero statistical aver¬ 
age. Notice that since the density is positive everywhere, 
nc{r) > 0, the field r]{r) is bounded above, ri{r) < /x-f A. 
We then insert Eq. (2) into the Gross-Pitaevskii equa¬ 
tion ( 1 ), and linearize it, which yields 

,(t) = (M + A) (3) 


perturbation theory [20, 21]. For stronger disorder, A is 
finite, and our approach accounts for the local depletion 
of the density around the disorder maxima. 

Knowing the density field ?Xc(r), we now treat the 
collective excitations. The phase and density fluctua¬ 
tions, 9 and 6n, are readily found by developing the 
many-body Hamiltonian up to order two in the opera¬ 
tor B{r) = (5ft(r)/2-y/nc(r) +‘i\/nc{r)9{r). The resulting 
quadratic Hamiltonian is then diagonalized by the Bo- 
goliubov transform B{r) = + where 

be is the annihilation operator of an elementary pair ex¬ 
citation of energy e. It yields the Bogoliubov-de Gennes 
equations [19, 22] 



where 


where, in Fourier space. 


i>(q) 


^(q) 

i + aiqP' 


( 4 ) 


The quantity F(r) is a generalized smoothed poten¬ 
tial [ 20 ] , where the healing length is renormalized by the 
shift A to the value ■\/ 477 i(/x -|- 3A/2). The zero- 

average condition, {r]{r)) = 0 , then yields 


0 = A-1-(min{F(r), 7 x-I-A/2}), (5) 

where (...) denotes statistical averaging. Note that 
Eq. (5) ensures that /x -|- 3A/2 > 0, so that is well 
dehned. In general, the density field nc(r) is therefore 
found by solving self-consistently Eqs. (4) and (5) for A 
and V{r), and Eq. (3) for ? 7 (r). 

This self-consistent solution is in good agreement with 
the exact numerical solution of the GPE (I) (see Fig. I). 
As expected, the density modulations r]{r) follow those 
of the disorder, smoothed at the length scale of the 
healing length [20, 21]. However, both the amplitude 
of the smoothed potential and the healing length are 
renormalized by the energy shift A. Moreover, for 
strong disorder, the field ? 7 (r) locally saturates to the 
constant value /x -|- A at positions where V (r) typically 
exceeds the chemical potential /x (more precisely, where 
14(r) > /x-|-A/2). In those regions, which will be referred 
to as depleted regions, we thus have nc(r) Ri 0 (see Fig. I). 
In order to interpret the shift A, we may rewrite Eq. (5) 
in the form A = dV PiV) \y — {p. + A/2)], where 

P{V) is the probability distribution of the smoothed 
potential and the integral is restricted to the depleted 
regions. The quantity A can thus be assimilated to the 
weight of the part of the smoothed potential that is 
truncated in the depleted regions. In particular, in the 
case of weak disorder for which V (r) never exceeds /x, we 
find A = 0 and we recover the solution given by usual 
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and 


= /^+^W-2r;(r) 

^ ' y +X7(r) — P(r)-I-2?7(r] 
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In this form, Eq. ( 6 ) devises a well-defined two-wave 
scattering problem. The dynamics of a given excitation 
at energy e is governed by the homogeneous propagator 
An and scattering from a disordered medium defined by 
U{y). The latter combines the two random fields V{r) 
and 7 ?(r), which are strongly correlated [see Eqs. (3) to 

( 5 )]. 

At this point, one could wonder whether the Bogoli- 
ubov approach is valid even in the presence of strong den¬ 
sity depletion. The main approximation here is the trun¬ 
cation of the many-body Hamiltonian at second order in 
the Bogoliubov operator. Since the latter is equivalent to 
the linearization of the time-dependent Gross-Pitaevskii 
equation (tGPE) [23] , it can be tested by comparing the 
meanfield dynamics predicted by the exact tGPE on the 
one hand and by the linearized tGPE on the other hand. 
Our results show excellent agreement between the two 
in all regimes, namely the phonon, particle, intermediate 
regimes for weak to strong density depletion (for details, 
see appendix. A). It validates the use of the Bogoliubov 
approach used here. 


B. One-parameter scaling theory 

Universal transport properties can now be inferred 
using the one-parameter scaling (OPS) approach [24], 
which can be extended to the case of excitations, as we 
outline here. It consists in developing a renormalization- 
group (RG) analysis of the size-dependent conductance. 
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The latter is identified to the Thouless number [25], 
which is the ratio of the energy scale associated to dif¬ 
fusion across a finite sample of size L, 6e = hDn/L'^ 
(with = Wels/d the classical diffusion constant, = 
h~"^\de/d\i\ the excitation velocity, and Zb the Boltzmann 
transport mean free path), to the energy-level spacing, 
Ae = \/N{e)L'^ (with N{e) the density of states per unit 
volume). In diffusive regimes, if kg is the momentum as¬ 
sociated to the energy e, then N{e) oc k'^~^/\de/d'k\ = 
k^~^/hwe, so that G{L) oc {kgl^){kgLY~‘^, and (3 = 
d\ogG/dlogL ^ d — 2. In localized regimes, the con¬ 
ductance is exponentially small, G{L) ^ exp(—L/Li^e) 
with Lioe the localization length, and /3 ~ logG. For 
d < 2, (3{G) is strictly negative and G{L) always flows 
down to the localized regime. Then, all states are local¬ 
ized, with the localization length cx Zb in ID and 
log(Lioc/ZB) oc kglg in 2D. Conversely, for d > 2, the 
RG flow has an unstable fixed point at kgl^ ~ 1, known 
as the mobility edge or the Anderson localization transi¬ 
tion [24]. Since the above scaling laws are independent 
of the dispersion relation, these features are all universal, 
except the transition point, which is determined by the 
value of the inverse disorder parameter (IDP) kgl^. 


C. Disorder parameter 

In order to estimate the IDP for the scattering prob¬ 
lem ( 6 ), we follow the approach of Ref. [7, 10] and 
extend it to strong disorder with possible local den¬ 
sity depletion (see details in appendix B). In brief, we 
note that the homogeneous propagator Cq does not sup¬ 
port only plane-wave modes with momentum kg such 
that h^k^/2m = Y+ A)^ — {ji + 2A), but also 
evanescent modes of penetration length 7 “^ such that 
k?Yll2.m = -I- (/X -I- A)2 + (m + 2A). The latter en¬ 

sures that for a scattering length larger than the pene¬ 
tration length, the excitation modes {ug,Vg) can be de¬ 
composed into two fields (g+,g“), where the second one 
is enslaved by the first one. Retaining only the leading 
disorder terms, the behavior of the excitation is then en¬ 
tirely determined by the field gf , which fulfills the closed 
equation 

(r) = (r) + Vg{r)g+ (r) (7) 

where 

Ve(r) = P(r) -/(£:) 77 (r), ( 8 ) 

with /(e) = ^ so-called screened 

potential Ve(r) results from the competition of the bare 
disorder V (r) and the meanfield repulsive interaction, de¬ 
termined by the field g(r). This competition is strongly 
energy dependent due to the factor /(e). Equation (7) 
describes an equivalent scattering problem, which can 



Figure 2. Schematic view of the two-impurity model defined 
by Eq. (10). It is made of two types of impurities with, re¬ 
spectively, positive (-t-V(/) and negative (— V),”) amplitudes. 
The impurities are randomly and independently spread over 
in space. Each impurity is Gaussian shaped with a width ctr. 


now be solved by standard quantum transport the¬ 
ory [26]. In the on-shell approximation [27], it yields 


1 2Trm? 

kglsie) ~ 


(1 - cos e)Cg [2kg sin( 6 l/ 2 )], 


( 9 ) 

where d^ld denotes the infinitesimal solid angle in d di¬ 
mensions and Ce(q) oc (|Ve(q)|^) is the power spectrum 
of the screened potential. Notice that in ID, the angular 
integral in Eq. (9) reduces to 0 = tt so that one recovers 
the result of Ref. [7, 10] for the Lyapunov exponent. 

It is worth pointing out that the previous ap¬ 
proach describes only excitations of energy e > Sg = 
a/2A(p. + 3A/2). Otherwise, we have < 0 and all 
modes of Co are evanescent. In the following, we will 
thus disregard the case of excitations £< Eg, which in 
most cases reduces to a very narrow energy range at the 
bottom of the spectrum, since A <C /r. It should as well 
be pointed out that since it retains only leading terms in 
disorder, our theory is not expected to be quantitatively 
exact, but rather to provide a qualitative description of 
the relevant physics. Possible extensions of the approach 
are discussed in the conclusion. 


III. LOCALIZATION OF BOGOLIUBOV 
QUASIPARTICLES IN AN IMPURITY MODEL 

A. The impurity model 

We can now discuss the behavior of the disorder pa¬ 
rameter of the excitations. For concreteness, let us con¬ 
sider a generic impurity model described by the potential 

3 

The impurities are independent Gaussian-shaped poten¬ 
tials of width (Tr, h{r) = exp(—r^/ 2 cr|), randomly dis¬ 
tributed in space with density p, and with amplitudes 
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Vj that take the values +Vf^ > 0 or —kg” < 0 with 
equal probability (see Fig. 2). The constant term V = 
— Vq)I2 ensures that the potential V{r) 
is of vanishing statistical average. The square disorder 
amplitude is = {V^) = p(v^CTr)'^[(Fo^) 2 + (yo“)^]/2. 
This two-impurity model generalizes the one-impurity 
model, which is widely used in studies of Anderson lo¬ 
calization in non-interacting systems [28, 29]. Here, 
we introduce two types of impurities, namely repulsive 
{Vj = > 0) and attractive {Vj = —V(^ < 0) ones. 

In contrast to non-interacting systems, it is crucial to dis¬ 
tinguish repulsive and attractive impurities in the present 
work because they have radically different effects on the 
density background. For instance, only repulsive impuri¬ 
ties can induce local density depletion. The two-impurity 
model is generic in the sense that it is the simplest one 
to describe a disordered system where different kinds of 
impurities are present. In addition, controlled impurity 
models can be realized in ultracold-atom systems where 
the impurities are made of individual atoms trapped at 
some random sites of an optical lattice [30]. So far, only 
one-impurity models have been realized [31] but they can 
be extended to models with different kinds of impurities 
using different atomic species. 


B. Behavior of the disordered parameter 

To compute the IDP for the previous impurity model, 
we first numerically determine the density background 
r]{r) and A following the previous self-consistent proce¬ 
dure. This permits to compute the screened potential 
Eq. (8) and its power spectrum, from which the IDP is 
inferred using Eq. (9). Figure 3 shows the energy de¬ 
pendence of the IDP for the 3D balanced impurity case 
(Pq^ = Vo“), plotted as a function of ke^- Qualitatively 
similar curves are found for lower dimensions and for im¬ 
balanced impurity cases (Pq^ ^ Depending on the 

disorder strength, the IDP exhibits three generic behav¬ 
iors. Notice that —>• 0 corresponds to e — 

For weak disorder (case A in Fig. 3), the IDP shows a 
non-monotonic energy dependence, which can be under¬ 
stood as follows. At high energy, the excitations are in¬ 
sensitive to the density background and behave as parti¬ 
cles in the bare disorder potential. Conversely, at low en¬ 
ergy, the excitations are strongly affected by the density 
background, which screens the disorder and suppresses 
scattering. More precisely, this holds when the chemical 
potential exceeds the maximum of the smoothed poten¬ 
tial, i.e. for V^h{0) — V<ji where h{r) is the smoothed 
impurity. Then, the density background as no depleted 
region. The power spectrum of the screened potential, 
Ce(q), can be computed explicitly as a function of that 
of the bare disorder, C'(q), and of the excitation energy 
e using Eq. (8). It yields 


Cs(q) 


A l + 4fcge^ 1 

^ 1 -k 2fc|^2 I q2^2 


C(q). 


( 11 ) 



Figure 3. Inverse disorder parameter (IDP) versus pair- 
excitation momentum for the 3D balanced impurity model, 
with (Tr = C pc-R = 2 X 10”“*, and for various values of the 
disorder amplitude Pq = Pq^ = Pq“. Shown are the results of 
Eq. (9), where the screened power spectrum Ce(q) is calcu¬ 
lated with the full screened disorder of Eq. (8) (dots) or with 
Eq. (11) (solid line), as well as the contributions of the bulk 
(dashed lines) and depleted (dotted lines) regions. 


Inserting this Eq. (II) into Eq. (9), we find the solid line 
in Fig. 3, which reproduces very well the data in the full 
energy range for case A. Notice that the same formula as 
found from Eqs. (9) and (11) was inferred from calcula¬ 
tions of the scattering mean free path using a different 
approach in Ref. [32]. Equation (11) generalizes the ID 
case [7, 10]. It defines a screening function [prefactor 
in the rhs of Eq. (11)], which can also be identified in 
single-scattering processes [21] and renormalizes the dis¬ 
order by the interactions. The behavior of the IDP can 
now be found by inspection of Eqs. (9) and (11). For 
kg ^ the screening is irrelevant and we can 

replace C£(q) by C'(q) in Eq. (9). We then recover the 
free-particle behavior [27, 33], 

kek (fceO^ kg » tX-'- (12) 

Conversely, for kg <C the screening strongly 

enhances kgl^^ compared to the free-particle case and we 
find 

(fceC)-", kg « r\^R'. (13) 

This result is in agreement with the universal behav¬ 
ior expected for Goldstone modes [18]. Both low-energy 
and high-energy scalings reproduce the behavior of case 
A and locate the minimum of the IDP at A:™'” ^ 
min(l/^, 1/(7r). Note that a non-monotonic behavior of 
the IDP is also found in the propagation of other kinds 
of waves, such as photonic [34] and acoustic [35] ones. 

For intermediate to strong disorder (cases B and C in 
Fig. 3), the energy dependence of the IDP found from the 
solution of the full scattering problem [Eqs. (8) and (9)] 
strongly differs from the weak disorder case at low energy, 
where kgl^ now increases with the energy. To understand 
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this, it should be noticed that, for h{0) — V > fi, the 
background density is now locally depleted around the 
positive impurities. Hence, during its propagation, an 
excitation goes through two types of regions, namely the 
density depleted region, and the rest, which constitutes 
the density bulk. In the bulk, the field ri(r) may be ap¬ 
proximated by H(r), provided we neglect the quantity 
A, which is valid for low impurity density, pa^ <g; 1. 
It yields a non-monotonic contribution to k^l^ similar 
to case A, with a smaller overall magnitude due to the 
truncation around the positive impurities (dashed lines 
in Fig. 3). Conversely, in the depleted regions, the field 
r]{r) saturates to the value p, + A. The bare disorder 
in those regions is thus protected against screening and 
Eq. (8) may be replaced by V£(r) ~ H(r) — (/i -f A)f{e). 
In this field, the excitations behave as (non-Goldstone) 
free particles, yielding a monotonic contribution to k^ls 
(dotted lines). In the white-noise limit {k^a-R <C 1), this 
contribution is 

(14) 

The various behaviors of the IDP observed in Fig. 3 can 
then be interpreted as follows. Neglecting the correla¬ 
tions between the contributions of the bulk and depleted 
regions, the disorder parameter (k^l^)~^ is approximately 
the sum of these two contributions. Its inverse (the IDP 
kjji, which is plotted in Fig. 3), is thus dominated by 
the smallest corresponding contribution. At low energy, 
because of the screening in the bulk, the contribution 
of the depleted region always dominates if it exists, and 
captures the free-particle-like behavior of k^l^. At inter¬ 
mediate energy the behavior of kels crucially depends on 
the relative magnitude of the two contributions. When 
V^h{0) — V > p (case B), only the upper fraction of 
the positive impurities is truncated and the bulk starts 
to dominate at moderate energy. It results in a turn¬ 
ing point where bulk and depleted regions equally 
contribute to the IDP, yielding there a local maximum. 
When increases, the depleted region dominates in 
a wider energy range and moves to higher values. 

When h{0) — V ^ p (case C), the positive impuri¬ 
ties are almost entirely truncated, and eventually 
merges with the local minimum which does not 

significantly depend on . The curve then becomes 
monotonic. 


C. Localization diagram 

We now turn to the localization properties of the collec¬ 
tive excitations, focusing on the 3D case where mobility 
edges appear at the localization threshold kj^ ^ 1- We 
determine the latter from IDP curves as those of Fig. 3 for 
the general model with different amplitudes of the posi¬ 
tive (Pq^) and negative (Pq”) impurities. In all cases, they 
are of type A, B, or C with an overall magnitude and posi¬ 
tions of the local maximum and minimum that depend on 
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Figure 4. Localization diagram of pair excitations in the 3D 
impurity model, plotted as a function of the amplitudes of 
positive (Vq"^) and negative (Pq~) impurities, for (Tr = ^ and 
po-R = 2 X 10~^. It exhibits four classes of mobility spec¬ 
tra, characterized by zero (0), one (I), two (II), or three (III) 
mobility edges. Note the different scales on the two axes. 



the parameters of the disorder potential. The resulting 
localization diagram, shown on Fig. 4, displays several 
localization classes with between zero and three mobil¬ 
ity edges in the excitation spectrum. The high energy 
states are always extended. The regions of the diagram 
are then determined by three conditions. Firstly, the ex¬ 
istence of depleted regions requires lg’^/i(0) — V > p, 
which defines the roughly vertical line on the diagram. 
On the left, the IDP curves are of type A and the low- 
energy excitations are extended {k^la > 1). Just on the 
right, they are of type B. The low-energy excitations are 
then localized {kj^ < 1), and there is at least one mo¬ 
bility edge. Secondly, when the local minimum of the 
IDP curve, (fc^B)™™, is below the localization thresh¬ 
old a band of localized states appear at intermediate en¬ 
ergy, giving rise to two additional mobility edges. For 
Vq <C Pq”, the condition reads 1 ^ (fcgZB)™^ oc 1 /(Pq“)^, 
which yields the nearly horizontal line on the diagram. 
Thirdly, when increases, the local maximum of the 
IDP curve decreases. In the region with three mobility 
edges (IB), the two low-energy ones disappear. In the 
region with one mobility edge (I), the IDP curve turn 
from type B to type C, without affecting the number of 
mobility edges. 

The localization diagram of Fig. 4 is expected to 
be generic. In particular, the competition of disorder, 
screening, and density depletion determine the diversity 
of mobility spectra. Yet, a given model of disorder does 
not necessary display all cases and the imbalanced impu¬ 
rity model with a finite correlation length is the simplest 
we found that does. For instance, for only positive im¬ 
purities or in the balanced case, the only possibilities 
are (0) or (I) because the minimum of the IDP can¬ 
not be controlled independently of the density depletion. 
Conversely, for only negative impurities, the depleted re¬ 
gion is absent and the only possibilities are (0) or (B). 
The case of white-noise disorder is also limited because 
the smoothed impurity potential diverges in the center, 
h{v) = e“’’/^/47r^^r, so that depleted regions strictly ap¬ 
pear as soon as Vq ^ 0, and the only possibilities left 
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are (I) and (III). 


IV. CONCLUSIONS 

The physics we have discussed here is particularly rel¬ 
evant to ultracold-atom experiments. In these systems, 
out-of-equilibrium dynamics can be generated by a local 
quench, which produces collective excitations [36]. In the 
presence of disorder, their transport properties and abil¬ 
ity to mediate long-range energy transfer are determined 
by the four classes of mobility spectra of the localization 
diagram. In case (0), all excitations are protected against 
localization and propagate diffusively, i.e. (r^) = 2D^t. 
In all other cases, energy can only be partially trans¬ 
ferred since some excitation modes are localized. Energy- 
resolved quenches may provide experimental evidence of 
such mobility spectra in ultracold gases. Moreover, these 
systems offer a wide range of models of disorder, e.g. im¬ 
purities [30, 31] and speckle potentials [2]. The statistical 
properties of the latter may be tailored, which may lead 
to even richer localization diagrams [33, 37, 38]. 

The observation of the localization effects we have 
discussed here requires that the lifetime r of the Bo- 
goliubov quasiparticles exceeds the transport meanfree 
time Tb = Ib/wc (with We the excitation group veloc¬ 
ity). On the one hand, for low temperature, the decay 
of Bogoliubov excitations is dominated by Beliaev pro¬ 
cesses [23, 39]. To estimate the corresponding decay rate 
r = l/r, we resort to local density approximation. The 
depleted regions, where the excitations behave as free 
particles with infinite lifetime, very weakly contribute 
to r. A good estimate of T is thus given by the bulk 
contribution. For a typical excitation e ~ it yields 


r ~ where = mg/4mh? is the scatter¬ 

ing length [39]. On the other hand, We — ^Jgn/m in the 
phonon regime and /b ~ 1/Ae in the region of interest, 
which yields Tb ~ m/A'Khnasc- Therefore, the validity 
condition of our approach reads Tb/t ~ ^/naZ ^ 1, 
which is the validity condition of the Bogoliubov ap¬ 
proach, very verified in dilute-gas Bose-Einstein conden¬ 
sates [23]. For instance, using the parameters of Ref. [40], 
we find find r ~ 6s and Tb ~ Sms, making experimental 
observation of our predictions possible. 

The approach used in this work provides a intuitive 
understanding of the physics at stake as well as generic 
qualitative predictions for the localization behavior of 
collective excitations. However, since it relies on lowest- 
order perturbation theory, it is not expected to be quan¬ 
titatively accurate. In particular, it does not take into 
account the disorder-induced shift of the dispersion re¬ 
lation, which is known, in the non-interacting case, to 
result in a possibly significant shift on the position of 
the mobility edge [33, 41]. Moreover, numerical calcu¬ 
lations in the non-interacting case show that the mo¬ 
bility edge significantly depends on the model of disor¬ 
der [42, 43]. Therefore, the determination of the precise 
localization diagram for collective pair excitations in dis¬ 
ordered Bose superfluids, as well as the identification of 
the various classes of mobility spectra predicted in the 
present work, require a full numerical resolution of the 
localization problem in the two-impurity model as well 
as in other models of disorder. 
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Appendix A: Validity of the second-order 
development in the presence of strong density 
depletion 

In this appendix, we show that the Bogoliubov ap¬ 
proach used in this paper, i.e. a development of the 
many-body Hamiltonian around the inhomogeneous den¬ 
sity background up to second order in fluctuations terms, 
is valid even in the presence of strong local density de¬ 
pletion. We perform this check by comparing exact cal¬ 
culations and linearized theory at the meanfield level. 
We consider the time-dependent Gross-Pitaevskii equa¬ 
tion (tGPE), 


ihdtif 


2m 


Ip -\-V{r) -i- g\tp\'^tp, 


(Al) 


which governs the time-evolution of a condensate wave- 
function ip(r,t). In the linearized approach, one writes 
ip{r,t) = y/njr) J- S'tp{r,t), where nc(r) is the density 
background found from the solution of the stationary 
GPE (1) and Stp{r,t) is a small perturbation. At low¬ 
est order, it yields the linearized equation 

(Pp.) = Ca. ( *f.) . (A2) 


where the matrix 

/ 

^GP = 

V 


2m 


V -b gric 
-gnc 


2m 


gric 


-V - gnc + pL) 


is exactly the one appearing in the Bogoliubov-de-Gennes 
equations [23]. This linearization procedure thus turns 
out to be equivalent to the Bogoliubov development of 
the many-body Hamiltonian to second order. Therefore, 
to check the validity of the latter in the presence of den¬ 
sity depletion, one can compare the results of the time- 
evolution of an excitation 5'ip{r,t) on top of a depleted 
condensate nc(r), using either the exact tGPE (Al) or 
its linearized version, Eq. (A2). 

We have performed this test in one dimension, which 
is the most unfavorable dimension due to large fluc¬ 
tuations. For the sake of simplicity, the density de¬ 
pletion is induced by a single strong potential barrier 
V{x) = Vb0(a — ja; — a;*]), of height Vq = 30g, width 
2a = 0.05L, and centered on some position a:*. The 
quantity 0 denotes the Heaviside function. In such a 
configuration, we first determine the density background 
nc(x) solving the stationary GPE Eq. (1) by imaginary 
time propagation. The latter is strongly depleted under 
the barrier (see Fig. 5). We then add a small excitation 
on top on this background, Sipo^x) = . 

The latter describes a plane wave of momentum k inside 
a Gaussian envelop of width a. In practice, we choose 
fc(T ^ 1 so that it is sufficiently narrow around k in mo¬ 
mentum space. We then compute the time evolution of 
this initial excitation either solving Eq. (Al) for tp{x,t) 







X 


-L 0 L 


x" 

3^ 


-L 0 L 

X 



X 

Figure 5. Time evolution of an excitation 5tpo{x) = 
gikx^-(x-XQ) jia on top of a depleted condensate in the pres¬ 
ence of a strong potential barrier V{x) = Vb0(a — \x — a:*|), 
for Vo = 30/r, 2a = 0.05L, and a = 0.05L. Three different 
regimes of the excitation spectrum are considered: (i) = 1 

with almost total reflection (top), (ii) = 5 with reflec¬ 
tion and transmission of the excitation of the same order to 
magnitude (central), (hi) = 10 with almost total trans¬ 
mission (bottom). Note that for clarity purposes, ^ is varied 
from a panel to another. Shown are the initial wavefunction 
‘ip{x, 0) = yjnc{x) -1- Sipoix) (solid red line) as well as the final 
wavefunction as given by the full time-dependent tGPE (solid 
orange line) and by the linearized tGPE (blue dotted line). 
The agreement between the two calculations is excellent in all 
cases. 


X 



kt-1 

t= 100, from (tGPE) 
g/L=0.005 t=100. from linearized (tGPE) 



with the initial condition 1 / 1 ( 0 ;, 0) = or 


solving Eq. (A2) for 6ip{x,t) with the initial condition 
6'ip{x,0) = 6'ipo{x). We have performed this compari¬ 
son in a wide range of parameters, from non-depleted to 
strongly depleted cases, and from the phonon <C 1) to 
the free-particle {k^ ^ 1) regimes. As shown on Fig. 5 in 
the case of strong depletion, we fonnd an excellent agree¬ 
ment in all cases, irrespective to the values of k and to 
the relative strength of reflection and transmission by the 
barrier. This validates the use of the linearized equation, 
and thus of the Bogoliubov approach, to study the dy¬ 
namics of the collective excitations even in the presence 
of strong modulations of the potential. 


Appendix B: Derivation of the inverse disorder 
parameter 


The background density field nc(r) being given by 
the solution of Eqs. (2) to (5), the collective excita¬ 
tions {ue-,Vs) can now be determined by solving the 
Bogoliubov-de Gennes equations (6), 


^0 





(Bl) 


where 


^0 — 


f/2m -b /i -b 2A 

V A 


-b/i “b A 
+h'^V‘^/2m — ^ 



and 


U{r) 


f+V{r) - 2'n{r) -77(r) \ 

-b77(r) -V^(r)-b 277(r)yl ■ 


The dynamics of a given excitation at energy e is thus 
governed by the homogeneous propagator £0 and scat¬ 
tering from the disordered medium defined by U (r). 

In order to solve the Bogoliubov-de Gennes equations 
(BdGEs), we generalize the approach of Ref. [7] to the 
strong disorder case where A yb 0. We first rewrite the 
BdGEs (Bl) in the form 



/—£ -b /r -b 2A /i -b A \ 
Y /x-bA e-b/x-b 2A J 



(B2) 


fV{r)-2r]{r) -r]{r) \ fuA 

V V' 


A suitable basis to perform diagrammatic expansion in 
leading disorder terms is found by applying the lin¬ 
ear transform that diagonalizes the 

homogeneous term in Eq. (B2), i.e. the matrix M = 
—£ -b /X “b 2A /X -b A 




A 


£ -b /X “b 2A 


. It yields 



/hWe 



— A-b£ 
-b A —£ 


2m 

2m 


-bA-b£ 
— A —£ 



,(B3) 
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where —h?k'^/2m = --y/eMT^^T+'Ap' + (/i + 2A) and 
= y^£ 2 + (/r + A)2 + (/r + 2A) are the eigen¬ 
values of the homogeneous matrix M. Without any ap¬ 
proximation at this stage, the BdGEs in the (g+, 5 “) 
basis then read 


2m 


W(r) = - 7 ;—V V{r) - f_{e)r]{r) g+{r) 


2m 

+<^+{e)g{r)g-{r) (B4) 

-srW =-Tj—V^5+(r)-h V{r) - f+{e)g{r) g-{r) 


2m 2m 

-h$-(e)? 7 (r) 5 +(r), 


(B5) 


with /±(e) = ^ ^ ^ ^ 

■y/e -K^A) ±(/j+A) ^ absence of disorder, Eqs. (B4) 

and (B5) are now decoupled and are straightforward to 
solve. The g'^ modes are plane waves of momentum k^, 
while the g~ are evanescent waves of penetration length 
7 “^. The latter vanish identically if the system is infinite 
or has periodic boundary conditions. 

In the presence of disorder, we can therefore make the 
assumption \g~\ <C \gf\ since g~ is at least one order of 
magnitude smaller that g'^ in and A/^. Keeping 

only the leading-order terms in Eq. (B5), we then find 


glected, which yields a closed equation for 

-^at (r) = (r) + Ve(r)g+ (r), (B7) 

where 


'^e{r) = V{r)-f{e)g{r) (B8) 

and, for simplicity, we now denote by /(e) the quantity 

/_(e) = quantity V£(r) defines 

a so-called screened potential of zero-average. It can be 
viewed as the screening of the bare potential V(r) by the 
density background encoded in g(r). It notably depends 
on the energy e of the Bogoliubov excitation. 

Therefore Eq. (B7), together with Eq. (B8), con¬ 
tains the leading disorder terms. It features an effec¬ 
tive single-wave scattering problem, which can now be 
solved by standard quantum transport theory [26]. Lo¬ 
calization properties are then determined in a two-step 
process [44, 45]. Firstly, the transport meanfree path is 
calculated in the semi-classical approach where interfer¬ 
ence of multiple-scattering paths is neglected. Within the 
on-shell approximation, which amounts to assimilate the 
spectral function to the disorder-free one, diagrammatic 
theory yields Eq. (9), 


2 ? 7 ? f 

9e - J dr' Gi/^^{r-r')r]{r')g+{r'), (B6) 

where Gi/^^ is the Green function associated to the dif¬ 
ferential operator —-I- I. In Fourier space, it reads 
Gi/^^(q) = (27r)“'^/^/[I -|- j(3^] and in real space it is 
a positive function of integral unity, decaying exponen¬ 
tially on a length scale I/ye. Then, since 2mjh'^'y'^ < g, 
and j4>_(e)j < 1 for any energy e, we find 

l5e"l^ J Gi/jei^-^')\v{^')\\9ti^')\dr'/g<{Vn/g)\g+\, 

which is consistent with the working assumption \g~\ <C 
jg/j. Therefore, the last term of Eq. (B4) can be ne¬ 


1 27rm^ 

kelsie) ~ 


(1 - cos e)Ce [2k, sin(6l/2)], 


(B9) 

for models of disorder with isotropic correlation func¬ 
tions [27], as considered in this work. For extension to 
anisotropic correlation functions, see Ref. [33]. Secondly, 
localization is found using either the one-parameter scal¬ 
ing theory [24] or the self-consistent approach [44, 45]. 
The one-parameter scaling theory is used and discussed 
in the main text. The self-consistent theory incorporates 
interference corrections on the top of diffusive dynamics, 
which yields a self-consistent equation for the diffusive 
constant or the localization length. Both approaches give 
the approximate localization threshold k,lB ^ 1 used in 
the present paper. 
























